R tiene una inmensa capacidad de graficar y visualizar datos de todo tipo, incluídos datos genéticos.
Las gráficas pueden hacerse desde la base de R (base) o con paquetes especializados en graficar, como lattice, o más recientemente ggplot2 y ggbio. También paquetes especializados en un tipo de datos que incluyen funciones para graficar, como ape para árboles filogenéticos.
En esta sección veremos una introducción a graficar en R usando graphics que es el sistema que viene con base y luego nos enfocaremos en gráficas más complejas y las principales usadas en análisis genéticos. En R se puede hacer mucho más que lo que veremos aquí, recomiendo profundizar.
Una de las mejores formas de aprender a hacer gráficas en R es buscar en internet/libro una gráfica parecida a la que queremos hacer y ver el código. Algunas recomendaciones:
Estas son las principales funciones para graficar utilizando la base de R. Puedes buscar ayuda de cada una con su nombre, y además en explorar argumentos extras con ?par
plot: generic x-y plottingbarplot: bar plotsboxplot: box-and-whisker plothist: histogramspie: pie chartsdotchart: cleveland dot plotsimage, heatmap, contour, persp: functions to generate image-like plotsqqnorm, qqline, qqplot: distribution comparison plotsDando x, y:
largo<-c(10,20,11,15,16,20)
ancho<-c(5,10,7,8,8,11)
plot(x=largo, y=ancho)
Dando un objeto que tiene dos columnas, se toman automático como x,y:
# ver el contenido de `cars`(una df ejemplo que viene con R)
head(cars)
## speed dist
## 1 4 2
## 2 4 10
## 3 7 4
## 4 7 22
## 5 8 16
## 6 9 10
plot(cars)
Si queremos especificar qué columnas serán x, y del objeto:
# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist)
Cambiar título de ejes e íconos:
# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist, xlab="Velocidad", ylab="Distancia", cex=0.5, pch=19)
Ejercicio: mira la ayuda de par y explica qué hacen los argumentos cex y pch.
Ejercicio: Repite la figura anterior pero cambiando los puntos por triángulos azules. Necesitarás esto.
Ejemplo con los datos islands (viene con R)
hist(islands)
Ejercicio En “Practicas/Uni7/data/reads.txt” encontrarás un archivo con la cantidad de lecturas de las muestras de tres librerías que fueron secuenciadas en Illumina. Grafica un histograma de las lecturas de cada muestra.
Ejemplo:
DNAcon<-data.frame(muestra=c("A", "B", "C"), concentracionADN=c(5,10,9))
barplot(DNAcon$concentracionADN, names.arg=DNAcon$muestra)
Ejercicio Repite la gráfica de DNAcon pero agregando títulos a los ejes x,y
Ahora cargemos un archivo real de datos:
reads<-read.delim("../Practicas/Uni7/data/reads.txt")
Hagamos una gráfica de barras y colorear acorde a info contenida en otra columna:
head(reads)
## Library sample nreads
## 1 Lib1 pobA21_r 1381230
## 2 Lib1 pobB10 1726622
## 3 Lib1 pobB05 819766
## 4 Lib1 pobC05 442508
## 5 Lib1 pobD21_r 1398874
## 6 Lib1 Outa112_r 2024684
barplot(reads$nreads, col=reads$Library)
Los colores que R ocupa para colorear algo están definidos en palette y pueden cambiarse
# Ver colores
palette()
## [1] "black" "red" "green3" "blue" "cyan" "magenta" "yellow"
## [8] "gray"
# Cambiar colores
palette(c("green", "blue", "red"))
# volver a graficar
barplot(reads$nreads, col=reads$Library)
Además de manualmente, los colores se pueden definir via paletas predeterminadas:
# Cambiar palette a 6 colores del arcoiris
palette(rainbow(6))
# volver a graficar
barplot(reads$nreads, col=reads$Library)
Checa otras palettes parecidas a rainbow en este link, y no te pierdas cómo nombrar muchos otros colores y utilizar otras paletas con más colores en la R Color Reference Sheet
# Graficar
barplot(reads$nreads, col=reads$Library)
# Agregar leyenda
legend(x="topleft", legend=levels(reads$Library), fill=palette()[1:3])
Nota que legend es una función por si misma (i.e. NO un argumento de plot) que requiere que antes de correrlo se haya corrido plot. Es decir una vez que creamos una gráfica podemos agregar sobre de esta una leyenda. Lo mismo puede hacerse con la función title.
Ejemplo:
boxplot(reads$nreads ~ reads$Library,
border = c("red", "blue", "darkgreen"))
Ejercicio agrega un título a la gráfica de boxplot de los reads.
La graficación de árboles filogenéticos se hace con el paquete ape, con el paquete phytools para funcionalidad extendida y con el paquete ggtree. Empezaremos por ape.
Bibliografía recomendada:
Los árboles filogenéticos pueden construirse en R, simularse en R o leerse a R.
Veamos un ejemplo con un árbol simulado:
# Cargar librería
library(ape)
## Warning: package 'ape' was built under R version 3.2.3
# Simular árbol
set.seed(1) # este comando es opcional, sirve para que todas "simulemos los mismos números" y podamos repetir de forma idéntica la simulación cada vez
tree <- rtree(n = 10, rooted=FALSE)
# ¿Qué tipo de objeto es?
class(tree)
## [1] "phylo"
# ¿Qué contiene?
tree
##
## Phylogenetic tree with 10 tips and 8 internal nodes.
##
## Tip labels:
## t10, t6, t9, t1, t2, t7, ...
##
## Unrooted; includes branch lengths.
str(tree)
## List of 4
## $ edge : int [1:17, 1:2] 11 12 13 13 12 11 14 15 15 14 ...
## $ tip.label : chr [1:10] "t10" "t6" "t9" "t1" ...
## $ edge.length: num [1:17] 0.718 0.992 0.38 0.777 0.935 ...
## $ Nnode : int 8
## - attr(*, "class")= chr "phylo"
## - attr(*, "order")= chr "cladewise"
# Graficar el árbol
plot.phylo(tree, edge.width=2)
La función plot.phylo puede abreviarse como plot. R sabe que debe usar la función plot.phylo y no plot básico (como lo usamos arriba) porque el objeto que le damos es un árbol.
Podemos modificar este árbol de manera similar a cómo lo hicimos en las gráficas anteriores.
Ejercicio Revisa la ayuda de plot.phylo y utiliza un argumento de dicha función para graficar el árbol simulado pero que se vea como un abanico y luego como un cladograma (así):
# plot árbol sin enraizar
plot.phylo(tree, edge.width=2)
# especificar output para enraizar:
tree<-root(tree, outgroup="t2")
# plot árbol enraizado
plot.phylo(tree, edge.width=2)
Ejercicio: Sigue los ejemplos de Phylogenetic trees in R, del blog Sensory Evolution para realizar los siguientes cambios al árbol anterior:
Podemos leer árboles en formato newick o nexus a R con la función read.tree de ape:
# cargar archivo
maiz.tree<-read.nexus("../Practicas/Uni7/data/tree")
# checar contenido
maiz.tree
##
## Phylogenetic tree with 165 tips and 163 internal nodes.
##
## Tip labels:
## maiz_3, maiz_68, maiz_91, maiz_39, maiz_12, maiz_41, ...
##
## Unrooted; includes branch lengths.
# graficar
plot(maiz.tree, type="unrooted", edge.width=0.1, cex=0.5)
Vamos a poner colores de acuerdo a la Categoría de Altitud en vez de nombres de muestras.
### Graficar por Categorías Altitud
# leer info extra de las muestras
fullmat<-read.delim("../Practicas/Uni7/meta/maizteocintle_SNP50k_meta_extended.txt")
# ¿Cuántos colores necesito?
numcolsneeded<-length(levels(fullmat$Categ.Altitud))
palette(rainbow(numcolsneeded))
# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)
# Agregar tip labels que correspondan a las categorías de altitud
tiplabels(pch=20, col=fullmat$Categ.Altitud)
# legend
legend(x= "bottomleft", legend=levels(fullmat$Categ.Altitud), pch=19, col=1:numcolsneeded, cex=1, , bty="n")
Ejercicio: Colorea el árbol anterior por Raza (sin incluir una leyenda porque son demasiadas)
Observa que las muestras 162-165 corresponden a teocintles (Zea m. mexicana y Z. m. parviglumis)
fullmat[162:165, 1:4]
## OrderColecta NSiembra Origen Raza
## 162 NA 96 Jardin Etnobotanico Oaxaca Zea m. parviglumis
## 163 NA 203 Luis Eguiarte Zea m. mexicana
## 164 NA 911 Luis Eguiarte 1_1 Zea m. parviglumis
## 165 NA 9107 Luis Eguarte 10_7 Zea m. mexicana
Si quisiéramos colorear todas las tips de maíz, pero no las teocintle podemos especificar esto así:
# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)
# Agregar tip labels sólo a los maíces
tiplabels(tip=c(1:161), pch=20, col="black")
Ejercicio Grafica el árbol de maíces de manera que los teocintles sean cuadrados negros y lo s maíces círculos verdes, así:
En R pueden visualizarse mapas de muchas maneras y de hecho hasta hacer análisis complejos con datos raster, como simulaciones de cambio climático y modelos de distribución potencial de las especies. Aquí sólo cubriremos brevemente cómo graficar un shapefile y agregar puntos.
Carguemos uno de los principales paquetes para manipular mapas:
library(maptools)
## Warning: package 'maptools' was built under R version 3.2.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.2.5
## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()
La función readShapePoly de dicho paquete nos permite leer un shapefile de polígonos (para puntos hay que usar otra función parecida ¿cómo crees que se llame?).
Por ejemplo vamos leer a R y graficar el shaphile de las regiones biogeográficas de México:
# leer shapefile
biogeo<-readShapePoly("../Practicas/Uni7/data/rbiog4mgw/rbiog4mgw.shp")
# plot
plot(biogeo)
# colorear por bioregioón
palette("default")
plot(biogeo, border="grey", col=biogeo$PROVINCIA)
## Cambiar colores default a algo más bonito
#¿Cuántos colores necesito?
levels(biogeo$PROVINCIA)
## [1] "Altiplano Norte (Chihuahuense)"
## [2] "Altiplano Sur (Zacatecano-Potosino)"
## [3] "Baja California"
## [4] "California"
## [5] "Costa del Pacifico"
## [6] "Del Cabo"
## [7] "Depresion del Balsas"
## [8] "Eje Volcanico"
## [9] "Golfo de Mexico"
## [10] "Los Altos de Chiapas"
## [11] "Oaxaca"
## [12] "Peten"
## [13] "Sierra Madre del Sur"
## [14] "Sierra Madre Occidental"
## [15] "Sierra Madre Oriental"
## [16] "Soconusco"
## [17] "Sonorense"
## [18] "Tamaulipeca"
## [19] "Yucatan"
# Generar paleta con colores de RColorBrewer
# ver opciones de colores
library(RColorBrewer)
display.brewer.all()
# generar paleta
palette(c(brewer.pal(9, "Set1"), brewer.pal(10, "Set3")))
# plot
plot(biogeo, border="grey", col=biogeo$PROVINCIA)
legend("bottomleft", legend=levels(biogeo$PROVINCIA), bty="n", cex=.4, fill=palette())
Ejercicio: cambia el color de las provincias “Tamaulipecas” y “Costa del Pacífico” a otro color.
Es muy común tener las coordenadas x,y de nuestros puntos de muestreo en un archivo junto con el resto de la info de nuestras muestras. Por ejemplo en el caso de la info de maíces que hemos estado utilizando:
fullmat[1:5,c("Latitud", "Longitud")]
## Latitud Longitud
## 1 28.68578 -107.58300
## 2 19.29083 -97.92944
## 3 28.53703 -108.92467
## 4 16.03225 -92.97972
## 5 20.09111 -101.11889
Agregar esta info a un mapa en la misma proyección y sistema de coordenadas puede hacerse con la función points:
# plot map
plot(biogeo, border="grey", col=biogeo$PROVINCIA, lwd=0.8)
# agregar puntos
points(fullmat$Longitud, fullmat$Latitud, pch=19, col="black", cex=0.4)
Ejercicio: Baja un mapa (nivel nacional) que te interese del GeoPortal de la CONABIO, ploteálo y agrega los puntos del muestro de maíz, utilizando una forma de punto diferente para los teocintles, y que los puntos estén coloreados por CategAltitud.
Hay una librería para bajar mapas (estáticos) de GoogleMaps: RgoogleMaps
Ejemplo:
# Cargar librería
library(RgoogleMaps)
# Establecer limites deseados
lat <- c(13,30) # definir limites en y
lon <- c(-112,-93) #definir limites en x
center = c(mean(lat), mean(lon)) #donde queremos centrar
zoom <- 4 #zoom deseado (1 lo más alejado)
# obtener mapa
mygooglemap <- GetMap(center=center, zoom=zoom, maptype= "terrain")
# plot
PlotOnStaticMap(mygooglemap)
png(), jpg(), pdf() según el formato en que queramos guardar.dev.offpng("../Practicas/Uni7/out/arbol.png")
plot(maiz.tree, type="unrooted", show.tip=FALSE)
tiplabels(pch=20, col="green3", tip=c(1:161), cex=.5)
tiplabels(pch=15, col="black", tip=c(162:165), cex=.5)
dev.off()
## quartz_off_screen
## 2
En R studio darle “Export” en el panel de la imagen y seleccionar un nombre de archivo y demás características.
Las gráficas que hemos visto hasta ahora pueden verse un poco feas de inicio y puede tomar un rato y mucho código arreglarlas a algo hermoso. ggplot2 es un paquete que ahorra este trabajo y que ha comenzado a ser ampliamente adoptado.
ggplot2 construye gráficas “definiendo sus componentes paso a paso”.
Para poder usar ggplot2 se requiere que la data.frame esté en formato largo, como vimos cuando revisamos la función gather. Además de esos apuntes puedes revisar esto si te quedan dudas.
Términos importantes:
Ojo: Mucho mejor que ver la ayuda de cada función es ver la Documentación online de ggplot2 y este R Graphics Cookbook
ggplot la función principal donde se especifican el set de datos y las variables a graficar.
geom_point()geom_bar()geom_density()geom_line()geom_area()geom_histogram()
aes los estéticos que pondremos: forma, transparencia (alpha), color, relleno, tipo de línea, etc.
scales para especificar si los datos se graficarán de forma continua, discreta, logarítmica.
themes para modificar los elementos de la gráfica no relacionados con los datos, como el tipo de letra y el color del fondo.
# Cargar ggplot2
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.2.4
# Examinar datos pre-cargados
head(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa
# graficar
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) + geom_point()
Pregunta: ¿Qué hace el símbolo +? Nota que el código anterior tmb puede escribirse así:
myplot<-ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width))
myplot + geom_point()
Los colores y formas se cambian en aes:
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width, color= Species, shape=Species)) + geom_point()
Ya sea en el aes de la función inicial o dentro de los geoms (Nota que el tamaño no es un aes, sino un argumento de geom_point)
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) +
geom_point(aes(color= Species, shape=Species), size=3)
Si queremos quitar el fondo gris:
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) +
geom_point(aes(color= Species, shape=Species), size=3) +
theme_bw()
Aveces queremos graficar en páneles separados la misma info para diferentes tratamientos o especies. Por ejemplo:
ggplot(data=iris, aes(x=Sepal.Length, y= Sepal.Width)) +
geom_point() +
facet_grid(Species ~ .)
Ejercicio Pon color por especie a la gráfica anterior:
Ejercicio Repite la gráfica anterior pero para que se vea así:
Ejercicio Repite la figura anterior pero cambiando los labels para que digan “Ancho de sépalo” y “Largo de sépalo”, respectivamente. Debe verse así:
También podemos agregar el resultado de un modelo matemático, como una regresión lineal:
En este tipo de gráficas la altura de las barras puede significar dos cosas:
stat="bin" en los argumentos de geom_bar.stat="identity" en los argumentos de geom_bar.Ejemplo:
# Cargar archivo
reads<-read.delim("../Practicas/Uni7/data/reads.txt")
# plot
p <- ggplot(data=reads, aes(x=sample, y=nreads, fill=Library)) +
geom_bar(stat="identity")
p
# Rotar nombres muestras
p <- p + theme(axis.text.x = element_text(angle = 45, hjust = 1, size=6))
p
La gráfica anterior no es igual a la que hicmos con barplot con los mismos datos ya que ggplot2 grafica en el orden de los levels, en este caso:
head(levels(reads$sample))
## [1] "NegC1" "NegC2" "NegC3" "Outa112" "Outa112_r" "Outa113"
Forma de solucionarlo:
# Cambiar orden de levels:
reads$sample<-factor(reads$sample, levels = reads$sample[order(1:length(reads$sample))])
head(levels(reads$sample))
## [1] "pobA21_r" "pobB10" "pobB05" "pobC05" "pobD21_r" "Outa112_r"
# Graficar
# plot
p <- ggplot(data=reads, aes(x=sample, y=nreads, fill=Library)) + geom_bar(stat="identity")
p + theme(axis.text.x = element_text(angle = 45, hjust = 1, size=6))
Siguiendo con los mismos datos anteriores:
# plot
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p
# quitar leyenda
p + guides(fill=FALSE)
Los árboles filogenéticos también pueden graficarse estilo ggplot, con el paquete de Biocounductor ggtree, que veremos más adelante junto con otros paquetes.
Al igual que en base, en ggplot es posible cambiar los colores manualmente o cambiando la paleta.
Recomiendo buscar más información y ejemplos en esta excelente guía Cookbook-R colores en ggplot2.
Ejemplos:
Cambiar colores manualmente
# Crear paleta:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# Usar paleta en gráfica:
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_manual(values=c("red", "blue", "green"))
Cambiar la paleta
# Crear paleta apta para daltónicos:
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# Usar paleta en gráfica:
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_manual(values=cbPalette)
Usar una paleta de ColorBrewer
p <- ggplot(data=reads, aes(x=Library, y=nreads, fill=Library)) + geom_boxplot()
p + scale_fill_brewer(palette="Set1")
Utilizar gradientes de colores
# Generar datos
set.seed(1)
df <- data.frame(xval=rnorm(50), yval=rnorm(50))
# Plot
ggplot(df, aes(x=xval, y=yval, colour=yval)) + geom_point()
# Cambiar gradiente
ggplot(df, aes(x=xval, y=yval, colour=yval)) + geom_point() +
scale_colour_gradientn(colours=rainbow(6))